int rand_a_b(int a, int b){
    return rand()%(b-a) + a;
}

using namespace std;

double myCRS(double *param_min, double *param_max, double ini_val[nbsobol][nbpara], double *best_point, double fun(double *)) {
    
    
	// ----------------------------
	// II) Controlled Random Search
	// ----------------------------
	// CRS2 - as in Kaelo and Ali (2006) "Some variants of the controlled random search algorithm for global optimization"

	int iteration_count = 0;
	bool STOP = false;
	
	double acceptance_threshold = 0.00001; // if the difference between the highest valued points and the lowest one is less than 1000, we stop.
	double EVALUATION[nbsobol]; // evaluation of the function at the N points
	int index_low, index_high, index;
	
	// for sample generation:
	int sample_size; // for random selection
	// int generated_id, id_ok; // if we want selection without replacement
    int SAMPLE[nbpara+1];
    #if randomnum == 1
        std::uniform_int_distribution<> unif_dist(0, (nbsobol-1));
    #endif
//    std::random_device rd; // do not use default_random_engine for random number generation, use mersene twister mt19937 (with random_device for seeding differently everytime)
//    std::mt19937 gen(rd()); //random generator
    std::mt19937 gen(123); // seed 123 (want always the same seed for reproducibility (relaunching the code to follow the process))
//    srand(123);
    
    bool belong_set;
	double new_point[nbpara];
	double new_coordinate;
	
	double new_evaluation;
	
	double POINTS[nbsobol][nbpara], POINTS2[nbpara];
    
	// Step 0: Initialization
	// ----------------------
    
    for(int w = 0; w < nbsobol; w++) {
        for(int d = 0; d < nbpara; d++) {
            POINTS[w][d] = ini_val[w][d];
        }
    }
    
    
	// Evaluate the function at the N generated points
    char evalwrite [200], evalwritesobol[200];
    FILE *writeeval, *writeeval2, *writesobol;
    
    sprintf(evalwrite,"INPUT/SMM_evaluation/evaluateSOBOL.out");
    strncpy(evalwrite, evalwrite, sizeof(evalwrite) - 1);
    evalwrite[sizeof(evalwrite) - 1] = 0;
    
    writeeval=fopen(evalwrite, "w");
    setbuf ( writeeval , NULL );
    
	for(int i = 0; i < nbsobol; i++) {
        for(int d = 0; d < nbpara; d++) {
            POINTS2[d] = POINTS[i][d];
        }
		EVALUATION[i] = fun(POINTS2);
  
        fprintf(writeeval,"%20.15f\n",EVALUATION[i]);
	}
    fclose(writeeval);
    
    //printf("\n \n EVALUATION DONE for %d SOBOL", nbsobol); getchar();
    // END OF SOBOL == 0 (ALL EVALUATION ARE WRITTENT) //
    

    
	// ----------
	// BEGIN LOOP
	// ----------

	while(STOP == false) {

		// Step 1: Find the best and worst points
		// --------------------------------------

		// Find the best (x_l) and worst points (x_h) in the set, where the best point has the lowest function value (f_l) and vice versa. 
		// Stop if f_h - f_l < epsilon

		index_low = 0; index_high = 0; // to initialize the first time

		for(int i = 0; i < nbsobol; i++) {
			if(EVALUATION[i] <= EVALUATION[index_low]) { index_low = i; } // EVALUATE[index_low] = f_l
			if(EVALUATION[i] >= EVALUATION[index_high]) { index_high = i; } // EVALUATE[index_low] = f_h
		}
  
        // write the new sobol //
        sprintf(evalwritesobol,"INPUT/SMM_para/POINTS_iter%d_high%d_low%d.out",iteration_count,index_high,index_low);
        strncpy(evalwritesobol, evalwritesobol, sizeof(evalwritesobol) - 1);
        evalwritesobol[sizeof(evalwritesobol) - 1] = 0;
        
        writesobol=fopen(evalwritesobol, "w");
        setbuf ( writesobol , NULL );
        
        for(int i = 0; i < nbsobol; i++) {
            for(int d = 0; d < nbpara; d++) {
               fprintf(writesobol,"%20.15f\n", POINTS[i][d]);
            }
        }
        
        fclose(writesobol);
    
    

		// Step 2: Stopping rule
		// ---------------------
		// put the stopping rule there because done after recomputing the index for high and low

		if((EVALUATION[index_high] - EVALUATION[index_low]) < acceptance_threshold) { STOP = true; }


		// Step 3: Generate a trial point
		// ------------------------------

		belong_set = false; // to initialize, say that the point does not belong to the set.

		while(belong_set == false) { // as long as the generated trial point does not belong to the set, do NOT evaluate the function

			// 1) Randomly pick n+1 points
			// ---------------------------
			// remark: we are using CRS2 method, so we will always include the point with the lowest value in this set (as the first point).

			// just pick n+1 index among the set 1:N (i.e. 0:(N-1) on C++)
			// <=> Generate n+1 different index (for CRS) (in reality only five index generated, the first one being the lowest index)

			sample_size = 1;
			SAMPLE[0] = index_low; // always include it as first point (CRS2)
			// Sample WITH replacement:
			while(sample_size < (nbpara + 1)) {
            #if randomnum == 1
				SAMPLE[sample_size] = unif_dist(gen);
            #endif
            #if randomnum == 2
                SAMPLE[sample_size] = rand_a_b(0,(nbsobol-1));
            #endif
				sample_size = sample_size + 1;
			}


			// 2) Generate new trial point
			// ---------------------------

			for(int j=0; j < nbpara; j++) { // param by param
				new_coordinate = 0;
				for(int i=0; i < (nbpara + 1); i++) { // point in sample by point in sample
					if(i != (nbpara)) {
						new_coordinate = new_coordinate + 2*POINTS[SAMPLE[i]][j]/(nbpara); //  the first points: 2 * mean (centroid)
						// remark: /!\ do this /nb_param at the end, if we do it first it bugs (int takes over double and make everything bug)
					} else {
						new_coordinate = new_coordinate - POINTS[SAMPLE[i]][j]; // for the last point we substract its value
					}
				}
				new_point[j] = new_coordinate;
			}

//            for(int j = 0; j < nbpara; j++) {
//                printf("%d %f %f %f", j, new_point[j], param_min[j], param_max[j]); getchar();
//            }
            
			// 3) Check if the point is within our Set limits
			// ----------------------------------------------

			belong_set = true;
			for(int i=0; i < nbpara; i++) {
				if(new_point[i] < param_min[i]) { belong_set = false; } // does not belong to set in this case
				if(new_point[i] > param_max[i]) { belong_set = false; }
			}
            
			// if does not belong to set, will just try some other point
		}
        
        
        
        
        
        
        
		// Step 4: Evaluate function at this new point
		// -------------------------------------------
		new_evaluation = fun(new_point);
        //printf("%f", new_evaluation);getchar();

        char evalNEWwrite [200];
        
        sprintf(evalNEWwrite,"INPUT/SMM_evaluation/evaluateNEW_%d.out",iteration_count);
        strncpy(evalNEWwrite, evalNEWwrite, sizeof(evalNEWwrite) - 1);
        evalNEWwrite[sizeof(evalNEWwrite) - 1] = 0;
        
        writeeval2=fopen(evalNEWwrite, "a");
        setbuf ( writeeval2 , NULL );
        
        fprintf(writeeval2,"%20.15f",new_evaluation);
        
        fclose(writeeval2);
    
    
    
		// Step 5: Replace if the evaluation is below the highest one
		// ----------------------------------------------------------

		if(new_evaluation < EVALUATION[index_high]) {

			for(int j = 0; j < nbpara; j++) { // Replace point
				POINTS[index_high][j] = new_point[j];
			}
			EVALUATION[index_high] = new_evaluation; // Replace evaluation
		}


		// update iteration count
		iteration_count = iteration_count + 1; // remark: consider that each function evaluation, even if > max, is a new iteration


		// and then go try a new point
	}
    
    
    for(int d = 0; d < nbpara; d++) {
        best_point[d] = POINTS[index_low][d];
    }
    
    return EVALUATION[index_low];
	
}
